home *** CD-ROM | disk | FTP | other *** search
/ Skunkware 5 / Skunkware 5.iso / src / X11 / lib / asincos.c next >
C/C++ Source or Header  |  1995-05-03  |  5KB  |  177 lines
  1. /*
  2.  * Copyright (c) 1985 Regents of the University of California.
  3.  * All rights reserved.
  4.  *
  5.  * Redistribution and use in source and binary forms are permitted
  6.  * provided that the above copyright notice and this paragraph are
  7.  * duplicated in all such forms and that any documentation,
  8.  * advertising materials, and other materials related to such
  9.  * distribution and use acknowledge that the software was developed
  10.  * by the University of California, Berkeley.  The name of the
  11.  * University may not be used to endorse or promote products derived
  12.  * from this software without specific prior written permission.
  13.  * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
  14.  * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
  15.  * WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
  16.  *
  17.  * All recipients should regard themselves as participants in an ongoing
  18.  * research project and hence should feel obligated to report their
  19.  * experiences (good or bad) with these elementary function codes, using
  20.  * the sendbug(8) program, to the authors.
  21.  */
  22.  
  23. #include "hdr.h"
  24. #include <errno.h>
  25.  
  26. /* ASIN(X)
  27.  * RETURNS ARC SINE OF X
  28.  * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits)
  29.  * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85.
  30.  *
  31.  * Required system supported functions:
  32.  *    copysign(x,y)
  33.  *    sqrt(x)
  34.  *
  35.  * Required kernel function:
  36.  *    atan2(y,x) 
  37.  *
  38.  * Method :                  
  39.  *    asin(x) = atan2(x,sqrt(1-x*x)); for better accuracy, 1-x*x is 
  40.  *          computed as follows
  41.  *            1-x*x                     if x <  0.5, 
  42.  *            2*(1-|x|)-(1-|x|)*(1-|x|) if x >= 0.5.
  43.  *
  44.  * Special cases:
  45.  *    if x is NaN, return x itself;
  46.  *    if |x|>1, return NaN.
  47.  *
  48.  * Accuracy:
  49.  * 1)  If atan2() uses machine PI, then
  50.  * 
  51.  *    asin(x) returns (PI/pi) * (the exact arc sine of x) nearly rounded;
  52.  *    and PI is the exact pi rounded to machine precision (see atan2 for
  53.  *      details):
  54.  *
  55.  *    in decimal:
  56.  *        pi = 3.141592653589793 23846264338327 ..... 
  57.  *    53 bits   PI = 3.141592653589793 115997963 ..... ,
  58.  *    56 bits   PI = 3.141592653589793 227020265 ..... ,  
  59.  *
  60.  *    in hexadecimal:
  61.  *        pi = 3.243F6A8885A308D313198A2E....
  62.  *    53 bits   PI = 3.243F6A8885A30  =  2 * 1.921FB54442D18    error=.276ulps
  63.  *    56 bits   PI = 3.243F6A8885A308 =  4 * .C90FDAA22168C2    error=.206ulps
  64.  *    
  65.  *    In a test run with more than 200,000 random arguments on a VAX, the 
  66.  *    maximum observed error in ulps (units in the last place) was
  67.  *    2.06 ulps.      (comparing against (PI/pi)*(exact asin(x)));
  68.  *
  69.  * 2)  If atan2() uses true pi, then
  70.  *
  71.  *    asin(x) returns the exact asin(x) with error below about 2 ulps.
  72.  *
  73.  *    In a test run with more than 1,024,000 random arguments on a VAX, the 
  74.  *    maximum observed error in ulps (units in the last place) was
  75.  *      1.99 ulps.
  76.  */
  77.  
  78. double asin(double x)
  79. {
  80.     double s,t,one=1.0;
  81.  
  82.     if(NAN(x)) return(x);
  83.  
  84.     s=copysign(x,one);
  85.  
  86.     if (s > one ) {
  87. #ifdef mmax
  88.       return (infnan(EDOM));
  89. #else
  90.       double zero = 0.0;
  91.       errno = EDOM;
  92.       return (zero/zero);
  93. #endif
  94.     }
  95.  
  96.     if(s <= 0.5)
  97.         return(atan2(x,sqrt(one-x*x)));
  98.     else 
  99.         { t=one-s; s=t+t; return(atan2(x,sqrt(s-t*t))); }
  100.  
  101. }
  102.  
  103. /* ACOS(X)
  104.  * RETURNS ARC COS OF X
  105.  * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits)
  106.  * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85.
  107.  *
  108.  * Required system supported functions:
  109.  *    copysign(x,y)
  110.  *    sqrt(x)
  111.  *
  112.  * Required kernel function:
  113.  *    atan2(y,x) 
  114.  *
  115.  * Method :                  
  116.  *                  ________
  117.  *                           / 1 - x
  118.  *    acos(x) = 2*atan2(  / -------- , 1 ) .
  119.  *                        \/   1 + x
  120.  *
  121.  * Special cases:
  122.  *    if x is NaN, return x itself;
  123.  *    if |x|>1, return NaN.
  124.  *
  125.  * Accuracy:
  126.  * 1)  If atan2() uses machine PI, then
  127.  * 
  128.  *    acos(x) returns (PI/pi) * (the exact arc cosine of x) nearly rounded;
  129.  *    and PI is the exact pi rounded to machine precision (see atan2 for
  130.  *      details):
  131.  *
  132.  *    in decimal:
  133.  *        pi = 3.141592653589793 23846264338327 ..... 
  134.  *    53 bits   PI = 3.141592653589793 115997963 ..... ,
  135.  *    56 bits   PI = 3.141592653589793 227020265 ..... ,  
  136.  *
  137.  *    in hexadecimal:
  138.  *        pi = 3.243F6A8885A308D313198A2E....
  139.  *    53 bits   PI = 3.243F6A8885A30  =  2 * 1.921FB54442D18    error=.276ulps
  140.  *    56 bits   PI = 3.243F6A8885A308 =  4 * .C90FDAA22168C2    error=.206ulps
  141.  *    
  142.  *    In a test run with more than 200,000 random arguments on a VAX, the 
  143.  *    maximum observed error in ulps (units in the last place) was
  144.  *    2.07 ulps.      (comparing against (PI/pi)*(exact acos(x)));
  145.  *
  146.  * 2)  If atan2() uses true pi, then
  147.  *
  148.  *    acos(x) returns the exact acos(x) with error below about 2 ulps.
  149.  *
  150.  *    In a test run with more than 1,024,000 random arguments on a VAX, the 
  151.  *    maximum observed error in ulps (units in the last place) was
  152.  *    2.15 ulps.
  153.  */
  154.  
  155. double acos(double x)
  156. {
  157.     double t,one=1.0;
  158.  
  159.     if(NAN(x)) return(x);
  160.  
  161.     if((x > 1.0) || (x < -1.0)) {
  162. #ifdef mmax
  163.       return infnan(EDOM);
  164. #else
  165.       double zero=0.0;
  166.       errno = EDOM;
  167.       return(zero/zero);
  168. #endif
  169.     }
  170.  
  171.     if( x != -1.0)
  172.         t=atan2(sqrt((one-x)/(one+x)),one);
  173.     else
  174.         t=atan2(one,0.0);    /* t = PI/2 */
  175.     return(t+t);
  176. }
  177.